%%这时涡旋三维结构合成function
load('WODTSMEC.mat',"TSMEC")
key = 'interp_pt_anomaly';
TSMEC.TSprofile.batchinterp();
TSMEC.TSprofile = TSMEC.TSprofile.profileLimit('len',30,'effictive_depth',1000);


ArgoTwoR = TSMEC.profileDetectBymultiR("InTwoR");
type = arrayfun(@(x) x.type, ArgoTwoR.ineddyArray);
AE_idx = find(type == 1);
AEArgo = ArgoTempSaltProfileClass.copy(ArgoTwoR,'index',AE_idx);
AEArgo.anomaly_measure_by_car09();%这是气候态法 
%上一步需要car09气候态数据如果没有可以尝试注释步骤（邻近法）
% Argoout = TSMEC.profileDetectBymultiR("outTwoR");
% AEArgo.anomaly_measure_by_background(Argoout);
AEArgo.anomly_check();
AEArgo = AEArgo.profileLimit('key',key);


AEArgo.data_clean_by_key(key);
X = arrayfun(@(x) x.deltaXY(1), AEArgo.ineddyArray);
Y = arrayfun(@(x) x.deltaXY(2), AEArgo.ineddyArray);
scatter(X,Y);
x = -2:0.2:2;
y = -2:0.2:2;

[xi,yi,zi,vq] = ArgoTempSaltProfileClass.scatterInterp(AEArgo.profile,  [X;Y], [x;y], key,'len',5,'lambda',50);
% vq = smooth3(vq);
depth = 10:10:1000';
EDC = EddyDrawClass();
EDC.DepthProfileDraw(vq, depth, 'clim',[-2 2]);
save("CE.mat","vq");


% R_surface = vq(:,:,10)';
% figure
% % R_surface = flipud(R_surface);
% %% 给合成涡旋设置背景场2°*2°
% [xx, yy] = meshgrid(x,y);
% [Fx, Fy] = gradient(R_surface);
% [dx,~] =  gradient(xx);
% [~,dy] =  gradient(yy);
% R = 6371000; % 地球半径，单位为 m    
% g = 9.8;
% dx = dx*111000;
% dy = dy*111000;
% f = gsw_f(15);
% u = -g./f.*Fy./dy;
% v = g./f.*Fx./dx;
% 
% quiver(xx,yy,u,v,1.5);




